---
title: "Main Results (IV)"
author: "Andrei Munteanu"
date: "11/01/2022"
output:
  html_document:
    fig_width: 8 
    fig_height: 8 
  # html_document:
    # toc: true 
    # toc_float: true
---
<style type="text/css">
.main-container {
  max-width: 1800px;
  margin-left: auto;
  margin-right: auto;
}
</style>

```{r options, include=FALSE, echo=FALSE, cache=FALSE}
# wd_write<-"G:/My Drive/Research/20190300 Romania Bac/Graphs/06 IV"
# #setwd(wd_write)
# knitr::opts_knit$set(root.dir=wd_write)
knitr::opts_knit$set(root.dir=getwd())
knitr::opts_chunk$set(cache.lazy=FALSE,cache=TRUE,warning=FALSE,echo=TRUE,comment='',prompt=T)
options(modelsummary_format_numeric_latex='plain')
knitr::opts_template$set(
  regular=list(include=TRUE,echo=TRUE,cache=FALSE,warning=FALSE,message=FALSE),
  regular_cache=list(include=TRUE,echo=TRUE,cache=TRUE,warning=FALSE,message=FALSE),
  invisible=list(include=FALSE,echo=FALSE,cache=FALSE,warning=FALSE,results=FALSE,message=FALSE),
  muted=list(include=TRUE,echo=TRUE,cache=FALSE,warning=FALSE,results=FALSE,message=FALSE),
  latex=list(include=TRUE,echo=TRUE,cache=FALSE))



data_regression<-as.data.frame(data_regression)
data_regression$n_hs_town_group<-data_regression$n_hs_town
data_regression$n_hs_town_group[data_regression$n_hs_town>=4 & data_regression$n_hs_town<=15]<-"4-15"
data_regression$n_hs_town_group[data_regression$n_hs_town>15]<-"16+"
data_regression$n_hs_town_group<-with(data_regression, reorder(n_hs_town_group, n_hs_town))

data_regression2<-data_regression %>%
 mutate(dec_town=relevel(dec_town,ref="6")) %>%
 filter(school_change==F)

data_regression2<-data_regression2 %>%
    mutate(n_students_bin=cut(n_students_town_yr,seq(from=0,to=100000,by=50))) 

```

## {.tabset}

### Tide
```{r, opts.label='regular'}
model_tide<-feols(grad_perc~entrance_perc+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
                         n_students_bin+as.factor(an)+specializare_bac2+scoala_de_provenienta+town,
                       cluster=~town,
                       data=data_regression2 %>% filter(school_change==F))
summary(model_tide)
```

```{r, opts.label='regular'}
model_tide2<-feols(grad_perc~entrance_perc+dec_town*n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
                         as.factor(an)+specializare_bac2+scoala_de_provenienta+town,
                       cluster=~town,
                       data=data_regression2 %>% filter(school_change==F))
summary(model_tide2)
```
### Track
Run regression:
```{r, opts.label='regular'}
model_track<-feols(grad_perc~entrance_perc+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
        as.factor(an)+specializare_bac2+scoala_de_provenienta+town+n_students_bin|
        class_mean~dec_town*n_hs_town_group,
      cluster=~town,
      data=data_regression2 %>% filter(school_change==F))
stage1_track<-summary(model_track,stage=1)
vcov_cluster_track<-stage1_track$cov.unscaled

```

Show first and second stages

```{r, opts.label='regular'}
summary(model_track)
stage1_track
```

```{r, opts.label='regular'}
model_track_ols<-feols(grad_perc~entrance_perc+class_mean+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
        as.factor(an)+specializare_bac2+scoala_de_provenienta+town,
      cluster=~town,
      data=data_regression2 %>% filter(school_change==F))


```

Show first and second stages

```{r,cache=FALSE}
summary(model_track_ols)
```
 
 
```{r, opts.label='regular'}
sizes<-unique(data_regression2[order(data_regression2$n_hs_town_group),]$n_hs_town_group)

quantile<-'dec_town'
n_quantiles<-10
deciles<-c("dec_town1","dec_town2","dec_town3","dec_town4","dec_town5","dec_town7","dec_town8","dec_town9","dec_town10")
deciles<-c("dec_town6",deciles[deciles %in% colnames(vcov_cluster_track)])


fe_stage_1_track<-data.frame()
for (i in 1:length(sizes)){
  size<-sizes[i]
  print("size")
  print(size)
  for (q in deciles) {
    print("q")
    print(q)
    q<-as.numeric(gsub("dec_town","",q))
    if (q!=6 & size!=1){
      town_name <-paste0('n_hs_town_group',size)  
      quart_name <-paste0(quantile,q)
      quart_town_name <-paste0(quantile,q,':n_hs_town_group',size)
      names<-c('entrance_perc',town_name,quart_name,quart_town_name)
      #names<-c(quart_name,quart_town_name)
      vcov_temp<-vcov_cluster_track[names,names]
      
      fe_town<-summary(stage1_track)$coefficients[town_name]
      fe_quant<-summary(stage1_track)$coefficients[quart_name]
      fe_town_quant<-summary(stage1_track)$coefficients[quart_town_name]
      
      fe<-ifelse(is.na(fe_town),0,fe_town)+ifelse(is.na(fe_quant),0,fe_quant)+ifelse(is.na(fe_town_quant),0,fe_town_quant)
      
      #fe<-ifelse(is.na(fe_quant),0,fe_quant)+ifelse(is.na(fe_town_quant),0,fe_town_quant)
    } else if (q==6 & size!=1){
      town_name <-paste0('n_hs_town_group',size)  
      names<-c('entrance_perc',town_name)
      vcov_temp<-vcov_cluster_track[names,names]
      #vcov_temp<-0
      
      fe_town<-summary(stage1_track)$coefficients[town_name]
      
      fe<-ifelse(is.na(fe_town),0,fe_town)
      #fe<-0
    } else if (q!=6 & size==1){
      quart_name <-paste0(quantile,q) 
      names<-c('entrance_perc',quart_name)
      vcov_temp<-vcov_cluster_track[names,names]
      
      fe_quant<-summary(stage1_track)$coefficients[quart_name]
      
      
      fe<-ifelse(is.na(fe_quant),0,fe_quant)
    } else if (q==6 & size==1){
      fe<-0
      se<-0
    }
    fe_grade<-summary(stage1_track)$coefficients['entrance_perc']*
      mean(data_regression2[data_regression2$dec_town==q & data_regression2$n_hs_town_group==size & !is.na(data_regression2$entrance_perc),]$entrance_perc)
    fe<-fe+fe_grade
    
    q<-q
    size<-size
    #se<-sqrt(sum(vcov_temp))
    #sqrt(vcov_cluster_track[town_name,town_name])+sqrt(vcov_cluster_track[quart_name,quart_name])+sqrt(vcov_cluster_track[quart_town_name,quart_town_name])+
    #            2*(vcov_cluster_track[town_name,quart_name])+2*(vcov_cluster_track[quart_name,quart_town_name])+2*(vcov_cluster_track[town_name,quart_town_name])
    if (q!=6 | size!=1){
      #fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se)
      se<-sqrt(t(as.matrix(c(fe_grade,rep(1,dim(vcov_temp)[1]-1))))%*%vcov_temp%*%(c(fe_grade,rep(1,dim(vcov_temp)[1]-1))))
      fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se) 
    } else {
      #fe<-data.frame(fe=0,decile=q,n_hs_towns=size,se=0)
      se<-fe_grade*summary(stage1_track)$se['entrance_perc']
      fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se) 
    }
    
    #print('dim')
    #print(dim(fe)[2])
    fe_stage_1_track<-bind_rows(fe_stage_1_track,fe)
  }
}
fe_stage_1_track$fe<-fe_stage_1_track$fe-fe_stage_1_track[fe_stage_1_track$decile==6 & fe_stage_1_track$n_hs_towns==1,]$fe
```

Track 1st Stage Graph:

```{r, opts.label='regular'}
lb<-min(fe_stage_1_track %>% mutate(lb=fe-se*qnorm(0.975)) %>% ungroup %>% dplyr::select(lb))
ub<-max(fe_stage_1_track %>% mutate(ub=fe+se*qnorm(0.975)) %>% ungroup %>% dplyr::select(ub))
graph_track_1<-ggplot(data=fe_stage_1_track,aes(x=decile,y=fe,color=n_hs_towns))+
  geom_errorbar(aes(ymin=fe-se*qnorm(0.975), ymax=fe+se*qnorm(0.975)),  width=0.5,size=0.8,position=position_dodge(0.5))+
  geom_point(size=1,position=position_dodge(0.5))+
   theme(strip.text.x = element_text(size=18),
        axis.text.x=element_text(angle = 0, vjust=0.5,size=14),
        axis.text.y = element_text(size=14),  
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        legend.title=element_text(size=18),
        legend.text=element_text(size=18))+
  #ylim(-0.1,0.3)+
  xlab("Own Admission Score Decile (Town-Cohort)")+
  ylab("Peer Admission Score Percentile (National-Cohort)")+
  scale_x_continuous(breaks=0:10)+
  scale_y_continuous(breaks=seq(from = -20, to = 60, by = 5),limits=c(lb,ub))+ 
  theme(legend.position = "none")
graph_track_1
```
Differences:

```{r, opts.label='regular'}
#Bottom decile difference in peer means (across diff town sizes)
fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns==1,]$fe-fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns=='16+',]$fe
#Top decile difference in peer means (across diff town sizes)
fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns==1,]$fe-fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns=='16+',]$fe
```





Track 2nd Stage Graph

```{r, opts.label='regular'}
name<-"fit_class_mean"
fe_stage_2_track<-fe_stage_1_track

fe_stage_2_track$se_b<-summary(model_track)$se[name]
fe_stage_2_track$se_mu<-fe_stage_2_track$se
fe_stage_2_track$b<-summary(model_track)$coefficients[name]
fe_stage_2_track$mu<-fe_stage_2_track$fe


names<-c("n_hs_town_group2","n_hs_town_group3","n_hs_town_group4-15","n_hs_town_group16+")
town_size_fe<-data.frame(n_hs_towns=factor(c(1,2,3,'4-15','16+'),levels=c(1,2,3,'4-15','16+')),
                          n_town=c(0,summary(model_track)$coefficients[names]),
                         se_town=c(0,summary(model_track)$se[names]),
                         covb=c(0,vcov(model_track)[names,"fit_class_mean"]))

fe_stage_2_track<-fe_stage_2_track %>% left_join(town_size_fe)
fe_stage_2_track<-fe_stage_2_track %>% mutate(se_mub=sqrt(se_b^2*mu^2+se_mu^2*b^2+se_b^2*mu^2))
fe_stage_2_track<-fe_stage_2_track %>% mutate(fe=b*mu+n_town)
```




Differences:

```{r, opts.label='regular'}

#Difference at graduation between 1st decile students
fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns==1,]$fe-
  fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns=='16+',]$fe
#Difference at graduation between 10th decile students
fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_2_track$n_hs_towns==1,]$fe-
  fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_2_track$n_hs_towns=='16+',]$fe

#Difference at graduation between top and bottom decile students in large towns
fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_1_track$n_hs_towns=='16+',]$fe-
  fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns=='16+',]$fe
#Difference at graduation between top and bottom decile students in small towns
fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_1_track$n_hs_towns==2,]$fe-
  fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns==2,]$fe
```

### School

Run regression:

```{r, opts.label='regular'}
#rm(model_track,stage1_track)
gc()
model_school<-feols(grad_perc~entrance_perc+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
        as.factor(an)+specializare_bac2+scoala_de_provenienta+town|
        school_mean~dec_town*n_hs_town_group,
      cluster=~judet_bac+town,
      data=data_regression2 %>% filter(school_change==F))
stage1_school<-summary(model_school,stage=1)
vcov_cluster_school<-stage1_school$cov.unscaled

```

Show first and second stages:

```{r, opts.label='regular'}
summary(model_school)
stage1_school
```

```{r, opts.label='regular'}
#rm(model_track,stage1_track)
gc()
model_school_ols<-feols(grad_perc~entrance_perc+school_mean+dec_town+n_hs_town_group+dec_town*n_students_town_yr+Unemployment_hs_bac*dec_town+Wages_hs_bac*dec_town+drop_hs_hs_bac*dec_town|
        as.factor(an)+specializare_bac2+scoala_de_provenienta+town,
      cluster=~judet_bac+town,
      data=data_regression2 %>% filter(school_change==F))


```

Show first and second stages:

```{r, opts.label='regular'}
summary(model_school_ols)
```

 
```{r, opts.label='regular'}
sizes<-unique(data_regression2[order(data_regression2$n_hs_town_group),]$n_hs_town_group)

quantile<-'dec_town'
n_quantiles<-10
deciles<-c("dec_town1","dec_town2","dec_town3","dec_town4","dec_town5","dec_town7","dec_town8","dec_town9","dec_town10")
deciles<-c("dec_town6",deciles[deciles %in% colnames(vcov_cluster_school)])


fe_stage_1_school<-data.frame()
for (i in 1:length(sizes)){
  size<-sizes[i]
  print("size")
  print(size)
  for (q in deciles) {
    print("q")
    print(q)
    q<-as.numeric(gsub("dec_town","",q))
    if (q!=6 & size!=1){
      town_name <-paste0('n_hs_town_group',size)  
      quart_name <-paste0(quantile,q)
      quart_town_name <-paste0(quantile,q,':n_hs_town_group',size)
      names<-c('entrance_perc',town_name,quart_name,quart_town_name)
      #names<-c(quart_name,quart_town_name)
      vcov_temp<-vcov_cluster_school[names,names]
      
      fe_town<-summary(stage1_school)$coefficients[town_name]
      fe_quant<-summary(stage1_school)$coefficients[quart_name]
      fe_town_quant<-summary(stage1_school)$coefficients[quart_town_name]
      
      fe<-ifelse(is.na(fe_town),0,fe_town)+ifelse(is.na(fe_quant),0,fe_quant)+ifelse(is.na(fe_town_quant),0,fe_town_quant)
      
      #fe<-ifelse(is.na(fe_quant),0,fe_quant)+ifelse(is.na(fe_town_quant),0,fe_town_quant)
    } else if (q==6 & size!=1){
      town_name <-paste0('n_hs_town_group',size)  
      names<-c('entrance_perc',town_name)
      vcov_temp<-vcov_cluster_school[names,names]
      #vcov_temp<-0
      
      fe_town<-summary(stage1_school)$coefficients[town_name]
      
      fe<-ifelse(is.na(fe_town),0,fe_town)
      #fe<-0
    } else if (q!=6 & size==1){
      quart_name <-paste0(quantile,q) 
      names<-c('entrance_perc',quart_name)
      vcov_temp<-vcov_cluster_school[names,names]
      
      fe_quant<-summary(stage1_school)$coefficients[quart_name]
      
      
      fe<-ifelse(is.na(fe_quant),0,fe_quant)
    } else if (q==6 & size==1){
      fe<-0
      se<-0
    }
    fe_grade<-summary(stage1_school)$coefficients['entrance_perc']*
      mean(data_regression2[data_regression2$dec_town==q & data_regression2$n_hs_town_group==size & !is.na(data_regression2$entrance_perc),]$entrance_perc)
    fe<-fe+fe_grade
    
    q<-q
    size<-size
    #se<-sqrt(sum(vcov_temp))
    #sqrt(vcov_cluster_school[town_name,town_name])+sqrt(vcov_cluster_school[quart_name,quart_name])+sqrt(vcov_cluster_school[quart_town_name,quart_town_name])+
    #            2*(vcov_cluster_school[town_name,quart_name])+2*(vcov_cluster_school[quart_name,quart_town_name])+2*(vcov_cluster_school[town_name,quart_town_name])
    if (q!=6 | size!=1){
      #fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se)
      se<-sqrt(t(as.matrix(c(fe_grade,rep(1,dim(vcov_temp)[1]-1))))%*%vcov_temp%*%(c(fe_grade,rep(1,dim(vcov_temp)[1]-1))))
      fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se) 
    } else {
      #fe<-data.frame(fe=0,decile=q,n_hs_towns=size,se=0)
      se<-fe_grade*summary(stage1_school)$se['entrance_perc']
      fe<-data.frame(fe=fe,decile=q,n_hs_towns=size,se=se) 
    }
    
    #print('dim')
    #print(dim(fe)[2])
    fe_stage_1_school<-bind_rows(fe_stage_1_school,fe)
  }
}
fe_stage_1_school$fe<-fe_stage_1_school$fe-fe_stage_1_school[fe_stage_1_school$decile==6 & fe_stage_1_school$n_hs_towns==1,]$fe
```


School 1st Stage Graph:

```{r, opts.label='regular'}
lb<-min(fe_stage_1_school %>% mutate(lb=fe-se*qnorm(0.975)) %>% ungroup %>% dplyr::select(lb))
ub<-max(fe_stage_1_school %>% mutate(ub=fe+se*qnorm(0.975)) %>% ungroup %>% dplyr::select(ub))
graph_school_1<-ggplot(data=fe_stage_1_school,aes(x=decile,y=fe,color=n_hs_towns))+
  geom_errorbar(aes(ymin=fe-se*qnorm(0.975), ymax=fe+se*qnorm(0.975)),  width=0.5,size=0.8,position=position_dodge(0.5))+
  geom_point(size=1,position=position_dodge(0.5))+
      theme(strip.text.x = element_text(size=18),
        axis.text.x=element_text(angle = 0, vjust=0.5,size=14),
        axis.text.y = element_text(size=14),  
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        legend.title=element_text(size=18),
        legend.text=element_text(size=18))+
  #ylim(-0.1,0.3)+
  xlab("Own Admission Score Decile (Town-Cohort)")+
  ylab("Peer Admission Score Percentile (National-Cohort)")+
  scale_x_continuous(breaks=0:10)+
  scale_y_continuous(breaks=seq(from = -30, to = 30, by = 5),limits=c(lb,ub))+ 
  theme(legend.position = "none")
graph_school_1
```

Differences:

```{r, opts.label='regular'}
#Bottom decile difference in peer means (across diff town sizes)
fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns==1,]$fe-fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns=='16+',]$fe
#Top decile difference in peer means (across diff town sizes)
fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns==1,]$fe-fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns=='16+',]$fe
```


School 2nd Stage Graph


```{r, opts.label='regular'}
name<-"fit_school_mean"
fe_stage_2_school<-fe_stage_1_school

fe_stage_2_school$se_b<-summary(model_school)$se[name]
fe_stage_2_school$se_mu<-fe_stage_2_school$se
fe_stage_2_school$b<-summary(model_school)$coefficients[name]
fe_stage_2_school$mu<-fe_stage_2_school$fe

names<-c("n_hs_town_group2","n_hs_town_group3","n_hs_town_group4-15","n_hs_town_group16+")
town_size_fe<-data.frame(n_hs_towns=factor(c(1,2,3,'4-15','16+'),levels=c(1,2,3,'4-15','16+')),
                          n_town=c(0,summary(model_school)$coefficients[names]),
                         se_town=c(0,summary(model_school)$se[names]),
                         covb=c(0,vcov(model_school)[names,"fit_school_mean"]))

fe_stage_2_school<-fe_stage_2_school %>% left_join(town_size_fe)
fe_stage_2_school<-fe_stage_2_school %>% mutate(se_mub=sqrt(se_b^2*mu^2+se_mu^2*b^2+se_b^2*mu^2))
fe_stage_2_school<-fe_stage_2_school %>% mutate(fe=b*mu+n_town)


```

Differences:

```{r, opts.label='regular'}

#Difference at graduation between 1st decile students
fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns==1,]$fe-
  fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns=='16+',]$fe
#Difference at graduation between 10th decile students
fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_2_school$n_hs_towns==1,]$fe-
  fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_2_school$n_hs_towns=='16+',]$fe

#Difference at graduation between top and bottom decile students in large towns
fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_1_school$n_hs_towns=='16+',]$fe-
  fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns=='16+',]$fe
#Difference at graduation between top and bottom decile students in small towns
fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_1_school$n_hs_towns==2,]$fe-
  fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns==2,]$fe
```

### Latex 2nd
```{r, opts.label='latex'}
f_big<-function(x) format(x, big.mark=",", scientific=FALSE, nsmall=1,digits=1)
f <- function(x) format(round(x/1000000,1) , nsmall = 1)
options(modelsummary_format_numeric_latex = "plain")

#print stage 2 table
setwd(paste0(wd_code,"/06_figures_03_04_table_04_IV"))
fileConn<-file("06_table_04.txt")
writeLines(print(
  modelsummary(list(model_school_ols,model_school,model_track_ols,model_track,model_tide2),
                            statistic = "std.error",
             estimate="{estimate}{stars}",
              stars=c('^{*}'=0.1,'^{**}'=0.05,'^{***}'=0.01),
               gof_map=list(list("raw" = "nobs", "clean" = "N", "fmt" = f),
                          list("raw" = "r.squared", "clean" = "R$^2$",fmt="%.2f")),
              metrics="R2",
              output="latex",
             fmt=2,
                escape=F)
), fileConn)
close(fileConn)




r2(model_school_ols)[2]
r2(model_school)[2]
r2(model_track_ols)[2]
r2(model_track)[2]
r2(model_tide)[2]
r2(model_tide2)[2]
```

### Latex 1st
```{r, opts.label='regular'}

summary(model_school,stage=1)
summary(model_track,stage=1)
```

```{r, opts.label='latex'}
options(modelsummary_format_numeric_latex = "plain")
f_big<-function(x) format(x, big.mark=",", scientific=FALSE, nsmall=1,digits=1)
f <- function(x) format(round(x/1000000,1) , nsmall = 1)

#print stage 1 table (appendix)
setwd(paste0(wd_code,"/06_figures_03_04_table_04_IV"))
fileConn<-file("06_table_04.txt")
writeLines(print(
  modelsummary(list(summary(model_school,stage=1) ,
                        summary(model_track,stage=1)),
                            statistic = "std.error",
             estimate="{estimate}{stars}",
              stars=c('^{*}'=0.1,'^{**}'=0.05,'^{***}'=0.01),
               gof_map=list(list("raw" = "nobs", "clean" = "N", "fmt" = f),
                          list("raw" = "r.squared", "clean" = "R$^2$",fmt="%.2f")),
              metrics="R2",
              output="latex",
                escape=F)

), fileConn)
close(fileConn)


r2(summary(model_school,stage=1))[2]
r2(summary(model_track,stage=1))[2]

```



### Graphs: 2 Pannels
```{r, opts.label='regular'}
fe_stage_2_track$Level<-'Track'
fe_stage_2_school$Level<-'School'
fe_stage_2_master<-rbind(fe_stage_2_track,fe_stage_2_school)

fe_stage_1_track$Level<-'Track'
fe_stage_1_school$Level<-'School'
fe_stage_1_master<-rbind(fe_stage_1_track,fe_stage_1_school)

lb<-min(fe_stage_1_master %>% mutate(lb=fe-se*qnorm(0.975)) %>% ungroup %>% dplyr::select(lb))
ub<-max(fe_stage_1_master %>% mutate(ub=fe+se*qnorm(0.975)) %>% ungroup %>% dplyr::select(ub))
graph_2_pannel_stage_1<-ggplot(data=fe_stage_1_master,aes(x=decile,y=fe,color=n_hs_towns))+
  geom_errorbar(aes(ymin=fe-se*qnorm(0.975), ymax=fe+se*qnorm(0.975)),  width=0.5,size=0.8,position=position_dodge(0.5))+
  geom_point(size=1,position=position_dodge(0.5))+
  facet_wrap(~Level,ncol = 2,scales='free_y')+
      theme(strip.text.x = element_text(size=18),
        axis.text.x=element_text(angle = 0, vjust=0.5,size=14),
        axis.text.y = element_text(size=14),  
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        legend.title=element_text(size=18),
        legend.text=element_text(size=18))+
  xlab("Own Admission Score Decile (Town-Cohort)")+
  ylab("Peer Admission Score Percentile (National-Cohort)")+
  scale_x_continuous(breaks=0:10)+
  scale_y_continuous(breaks=seq(from = -30, to = 35, by = 5),limits=c(lb,ub))+
  labs(color=str_wrap("Number of High Schools in Town",width=10))+
  theme(legend.title.align=0.5)

setwd(wd_code)
setwd(paste0(wd_code,"/06_figures_03_04_table_04_IV/other output"))
pdf("figure_iv_1st_stage_2_panel.pdf" , width = 12 , height = 6)
  graph_2_pannel_stage_1
dev.off()

setwd(wd_code)
setwd(paste0(wd_code,"/06_figures_03_04_table_04_IV"))
pdf("06_figure_03.pdf" , width = 12 , height = 6)
  graph_2_pannel_stage_1 + scale_color_grey()+ scale_fill_grey()
dev.off()

lb<-min(fe_stage_2_master %>% mutate(lb=fe-se*qnorm(0.975)) %>% ungroup %>% dplyr::select(lb))
ub<-max(fe_stage_2_master %>% mutate(ub=fe+se*qnorm(0.975)) %>% ungroup %>% dplyr::select(ub))
graph_2_pannel_stage_2<-ggplot(data=fe_stage_2_master,aes(x=decile,y=fe,color=n_hs_towns))+
  geom_errorbar(aes(ymin=fe-se*qnorm(0.975), ymax=fe+se*qnorm(0.975)),  width=0.5,size=0.8,position=position_dodge(0.5))+
  geom_point(size=1,position=position_dodge(0.5))+
    facet_wrap(~Level,ncol = 2,scales='free_y')+
         theme(strip.text.x = element_text(size=18),
        axis.text.x=element_text(angle = 0, vjust=0.5,size=14),
        axis.text.y = element_text(size=14),  
        axis.title.x = element_text(size=14),
        axis.title.y = element_text(size=14),
        legend.title=element_text(size=18),
        legend.text=element_text(size=18))+
  facet_wrap(~Level,ncol = 2,scales='free_y')+
  xlab("Own Admission Score Decile (Town-Cohort)")+
  ylab("Own Graduation Score Percentile (National-Cohort)")+
  scale_x_continuous(breaks=0:10)+
  scale_y_continuous(breaks=seq(from = -40, to = 40, by = 2),limits=c(lb,ub))+
  labs(color=str_wrap("Number of High Schools in Town",width=10))+
  theme(legend.title.align=0.5)

setwd(wd_code)
setwd(paste0(wd_code,"/06_figures_03_04_table_04_IV/other output/"))
pdf("figure_iv_2nd_stage_2_panel.pdf" , width = 12 , height = 6)
  graph_2_pannel_stage_2
dev.off()

setwd(wd_code)
setwd(paste0(wd_code,"/06_figures_03_04_table_04_IV/"))
pdf("06_figure_04.pdf" , width = 12 , height = 6)
  graph_2_pannel_stage_2 + scale_color_grey()+ scale_fill_grey()
dev.off()

  graph_2_pannel_stage_2


```

These are the main results!

```{r, opts.label='regular'}
graph_2_pannel_stage_1
graph_2_pannel_stage_2
```




### Differences in Achievement:
Widening of achievement gap, top vs bottom students, same town size (large town):
```{r, opts.label='regular'}
fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_2_school$n_hs_towns=='16+',]$fe- fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns=='16+',]$fe
fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_2_track$n_hs_towns=='16+',]$fe- fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns=='16+',]$fe
```

Widening of achievement gap, top vs bottom students, same town size (small town):
```{r, opts.label='regular'}
fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_2_school$n_hs_towns=='1',]$fe- fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns=='1',]$fe
fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_2_track$n_hs_towns=='1',]$fe- fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns=='1',]$fe
```

<!-- fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_2_school$n_hs_towns=='2',]$fe-fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns=='2',]$fe -->
<!-- fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_2_track$n_hs_towns=='2',]$fe-fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns=='2',]$fe -->


Widening of achievement gap, top urban vs top rural students:
```{r, opts.label='regular'}
fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_2_school$n_hs_towns=='16+',]$fe- fe_stage_2_school[fe_stage_2_school$decile==10 & fe_stage_2_school$n_hs_towns=='1',]$fe
fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_2_track$n_hs_towns=='16+',]$fe- fe_stage_2_track[fe_stage_2_track$decile==10 & fe_stage_2_track$n_hs_towns=='1',]$fe
```

Widening of achievement gap, bottom urban vs bottom rural students:
```{r, opts.label='regular'}
fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns=='16+',]$fe- fe_stage_2_school[fe_stage_2_school$decile==1 & fe_stage_2_school$n_hs_towns=='1',]$fe
fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns=='16+',]$fe- fe_stage_2_track[fe_stage_2_track$decile==1 & fe_stage_2_track$n_hs_towns=='1',]$fe
```

### Differences in Sorting:
Peers of low-ability students (large vs small town):
```{r, opts.label='regular'}
fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns=='1',]$fe- fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns=='16+',]$fe
fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns=='1',]$fe- fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns=='16+',]$fe
```

Peers of high-ability students (large vs small town):
```{r, opts.label='regular'}
fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns=='1',]$fe- fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns=='16+',]$fe
fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns=='1',]$fe- fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns=='16+',]$fe
```

Peers of large town students:
```{r, opts.label='regular'}
fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns=='16+',]$fe- fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns=='16+',]$fe
fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns=='16+',]$fe- fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns=='16+',]$fe
```
Peers of small town students:
```{r, opts.label='regular'}
fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns=='1',]$fe- fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns=='1',]$fe
fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns=='1',]$fe- fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns=='1',]$fe
```

Sorting school vs track
```{r, opts.label='regular'}
fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns=='1',]$fe- fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns=='1',]$fe
fe_stage_1_track[fe_stage_1_track$decile==10 & fe_stage_1_track$n_hs_towns=='16+',]$fe- fe_stage_1_school[fe_stage_1_school$decile==10 & fe_stage_1_school$n_hs_towns=='16+',]$fe

fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns=='1',]$fe- fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns=='1',]$fe
fe_stage_1_track[fe_stage_1_track$decile==1 & fe_stage_1_track$n_hs_towns=='16+',]$fe- fe_stage_1_school[fe_stage_1_school$decile==1 & fe_stage_1_school$n_hs_towns=='16+',]$fe
```














